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Summary 

This report describes a novel numerical procedure 
for the iterative solution of inviscid flow problems 
and demonstrates its utility for the calculation of 
steady subsonic and transonic flow fields. The 
method is more general than previously developed 
iterative methods in that no assumptions concerning 
the existence of either a velocity potential function or 
a stream function are required. 

Application of the surrogate-equation technique 
defined herein allows the formulation of stable, fully 
conservative, type-dependent finite difference 
equations for use in obtaining numerical solutions to 
systems of first-order partial differential equations, 
such as the steady-state Euler equations. Included 
among the results presented are steady, two- 
dimensional solutions to the Euler equations for both 
subsonic, rotational flow and supersonic flow and to 
the small-disturbance equations for transonic flow. 
A computational efficiency in excess of that obtained 
by means of the standard perturbation-potential 
approach is indicated for the small-disturbance 
equations. Possible improvements to, and extensions 
of, the method are discussed. 


Introduction 


Motivating Factors 

The present study is concerned with the numerical 
solution of steady inviscid flow problems. Many 
important physical situations encountered in modern 
engineering and applied science can be accurately 
modeled within the constraints of steady inviscid flow 
theory. Timely substantiation of this claim is provided 
by the generally good results presently being obtained 
from the use of such mathematical models for the 
design and analysis of transonic airfoils. 

Most work, both past and present, has, however, 
dealt with that subset of inviscid flows that are 
irrotational and hence for which a velocity potential 
function exists. Although many flows of interest can 
be successfully modeled within this additional 
restriction, many others, constituting in all 
probability a larger class, cannot. Certainly all those 
flows of practical interest where significant gradients 
of entropy or total enthalpy can occur require a more 
general model than one based on potential flow 


theory. Significant gradients are virtually certain to 
occur in many internal flows, in particular in those 
through modern turbomachinery, and can also occur 
in external flows, especially when the flow over a 
number of interacting components is considered. 

Of course the hyperbolic partial differential 
equations describing supersonic steady inviscid flow 
problems can be solved, for both potential and more 
general flow situations, by means of existing 
mathematical and numerical techniques. 
Consequently such flows are not the object of the 
present study. Rather it is the subsonic and transonic 
flow problems, described respectively by elliptic and 
mixed elliptic-hyperbolic equations, with which the 
research described herein deals. 

At this juncture the issue of computational 
efficiency makes its importance felt. Subsonic and 
transonic inviscid flow problems can be solved by 
computing a temporal asymptote to the unsteady 
equations of motion. However, the computation 
times can be quite long. In contrast to this approach 
the present research describes a method for the direct 
solution of the steady equations. By proceeding in 
such a fashion the resolution of the transient physical 
states between the initial state and the desired final 
state is avoided. Thus a means for the more efficient 
solution of subsonic and transonic steady inviscid 
flow problems is provided. As described 
subsequently in this report the method is based on the 
creation of a higher order system that serves as a 
surrogate for the first-order partial differential 
equations of inviscid flow theory. 

Literature Review 

When written in primitive variable form the 
systems of partial differential equations used to 
describe the steady motions of an inviscid fluid are of 
first order and of mixed elliptic-hyperbolic type. 
Common examples of such systems are the transonic 
small-disturbance equations, the potential flow 
equations, and the Euler equations. 

Because of the difficulties associated with both the 
formulation of robust finite difference analogs for 
such equations and the construction of stable 
iterative procedures for their numerical solution 
(refs. 1 and 2), these partial differential systems are 
not usually solved in their steady, primitive variable 
form. Rather, as is well known, the transonic small- 
disturbance equations and the potential flow 
equations are transformed into scalar, second-order. 


panial differential equations by introducing either a 
velocity potential function (refs. 3 and 4) or a stream 
function (refs. 5 and 6). The steady Euler equations, 
on the other hand, are replaced by their unsteady 
versions, for which a temporally asymptotic steady 
solution is sought, either in real time (refs. 7 and 8) or 
in pseudo time (refs. 9 to 11). 

Relatively few departures from these approaches 
are to be found in the literature. Steger and Lomax 
(ref. 12) developed an iterative procedure for solving 
a nonconservation form of the steady Euler 
equations for subcritical flow with small shear. 
Chattot (ref. 13) solved the transonic small- 
disturbance equations by differentiating them to 
obtain a second-order system, a special case of the 
approach discussed herein. He later adopted a 
variational formulation (refs. 14 and 15) and applied 
it to model problems representing the Euler 
equations. Ozer (ref. 16) developed a relaxation 
procedure for solving the equations of motion when 
they are reformulated in such a manner as to yield a 
second-order partial differential equation in the 
logarithm of the pressure, together with first-order 
equations for the remaining variables. 

The work of these authors notwithstanding, 
contemporary numerical simulations of steady 
inviscid flow generally resort to either relaxation 
solutions of steady second-order equations with 
derived dependent variables or time-asymptotic 
solutions of unsteady first-order systems. In the 
former case, generality is lost; in the latter case the 
computational efficiency can be quite low. 

Scope of Present Study 

As is readily apparent in the foregoing literature 
survey there are two general approaches to the 
numerical solution of steady inviscid flow problems. 
One approach involves the time-accurate solution of 
the complete, unsteady Euler equations of motion. 
Taken in their time-dependent form the governing 
Euler equations are of hyperbolic type and their 
numerical solution is a relatively straightforward 
matter. Hence, one may attempt to obtain the 
solution to a steady flow problem as the temporal 
asymptote of the solution to an unsteady flow 
problem. This approach has been successfully 
employed by several researchers. 

An alternative approach is to solve the Euler 
equations by a method that is not time accurate but 
that produces the desired steady-state result. Such 
methods are generally referred to as relaxation or 
iterative methods. There is substantial opinion and a 
considerable body of evidence that relaxation 
methods can provide a converged solution more 
quickly than can time-accurate methods and hence 


lead to more efficient use of computer resources. 
Because of the limited capacity of presently available 
computers and the complex nature of the phenomena 
under investigation, this issue of computational 
efficiency is of great importance. This is especially 
true if numerical solutions are to be used for design 
purposes, which typically require large numbers of 
cases to be computed. 

At present, steady-state solutions of the Euler 
equations for subsonic and transonic flow problems 
are found primarily by means of the time-accurate 
approach. Since this approach can in the words of 
Lomax and Steger (ref. 2) be “disastrously slow,” 
the development of efficient relaxation procedures 
could be quite beneficial. It is, however, no 
coincidence that with the few exceptions indicated 
previously little has been achieved toward creating a 
suitable relaxation procedure for the complete, 
steady Euler equations. The task of constructing 
stable iteration algorithms for equations such as the 
Euler equations, which involve only first derivatives 
of primitive variables, has received little attention "in 
the literature. The questions as to whether, and under 
which circumstances, such procedures exist have not 
yet been satisfactorily answered. 

The present work describes a new procedure, 
referred to as the surrogate-equation technique, that 
is designed to circumvent the difficulties associated 
with the nature of the steady Euler equations and to 
permit their solution by means of conventional 
iterative techniques. The numerical stability 
problems normally encountered when attempting to 
formulate finite difference equations for the steady- 
state, first-order Euler equations in regions where 
they exhibit elliptic behavior are avoided by 
introducing an alternative higher-order partial 
differential system for which proven numerical 
solution techniques are readily available. 

For clarity the essential ideas comprising the 
surrogate-equation technique are introduced within 
the context of a quite simple and well-understood 
first-order system of partial differential equations, 
namely, the Cauchy-Riemann system. More 
interesting subsonic and transonic flow problems are 
subsequently discussed and solved. 

By proceeding in this fashion a basis for the 
iterative solution of subsonic and transonic inviscid 
flow problems that lie beyond the restrictions of 
potential flow theory is developed. Although the 
efficient use of computer resources is a motivating 
force, no attempt is made here to proceed beyond the 
now standard, successive-line-relaxation solution 
procedure. The application of convergence 
acceleration techniques is reserved for future study. 

It is nevertheless interesting to note that when applied 
to the transonic small -disturbance equations, as 


discussed in the section Transonic Flow, the 
surrogate-equation technique leads to an algorithm 
that, on the basis of the computational 
experimentation reported herein, appears to exceed 
the computational efficiency of the standard 
Murman and Cole algorithm by several multiples. 


Cauchy-Riemann Problem 

To illustrate the essential aspects of the surrogate- 
equation technique, we examine its application to a 
first-order system of partial differential equations of 
minimal complexity. Consider the equations 
describing in two dimensions the flow of an 
incompressible and irrotational fluid. 


fact that most of the equations of mathematical 
physics are second-order partial differential 
equations means that a large variety of proven 
numerical methods exist and are at our disposal. 

Surrogate-Equation Formulation 

Consider now another approach to obtaining a 
solution to equations (1) and (2) that for convenience 
is referred to as the surrogate-equation technique. This 
technique is of general application and is not restricted 
to the class of irrotational, incompressible fluid flows 
that we are treating here only by way of an example. 

By defining the two-component vectors / and g 
such that 


Ux + Vy=0 
Vx — Uy=0 



where u and v represent the components of velocity 
in the x and y directions, respectively. As is well 
known, these equations are simply the Cauchy- 
Riemann equations. 


and 



Potential Formulation 

To solve equations (1) and (2), hydrodynamicists 
have long made use of the fact that a velocity 
potential, <p = <p{x,y), can be introduced such that the 
condition of irrotationality (eq. (2)) is identically 
satisfied: 

<Px = u 


we can rewrite equations (1) and (2) as 

fx+gy=0 (3) 

Equation (3) can then in turn be reexpressed as 
fx +Afy =0, where the Jacobian matrix A is defined 
such that A = dgldf. Since A is a. constant matrix, we 
can write fx + (.AJ)y = 0 or more conveniently in 
differential operator notation. 


iPy=V 

Substitution of the velocity potential into the 
incompressibility condition (eq. (1)) then leads to a 
Laplace equation in <p: 

'Pxx + ‘Pyy = 0 

A succinct and informative discussion of this 
development from the viewpoint of complex variable 
theory is given in reference 17. 

Hence for a particular incompressible and 
irrotational flow problem the resulting boundary 
value problem for the first-order Cauchy-Riemann 
system can be reformulated as a boundary value 
problem for the Laplace equation in the single scalar 
dependent variable <^. As the theory of harmonic 
functions is a quite mature branch of mathematics, a 
large number of particular solutions are available for 
this problem. Furthermore, if circumstances should 
indicate the desirability of a numerical solution, the 




(4) 


where I represents the two-dimensional identity 
matrix. 

If we now operate on equation (4) with the 
differential operator 


with M=I and N=—A (or N=A^), we obtain a 
second-order partial differential equation for / that 
has many pleasing properties. The form of M and N 
might be suggested by analogy with holomorphic 
function theory (ref. 18). In any case the present 
exposition undertakes to illustrate the utility of such 
a choice. 
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Clearly the equation resulting from an application 
of operator (3) to equation (4) is 

In the case under consideration, where 



and hence equation (6) reduces quite 

simply to 



This is the two-dimensional Laplace equation for the 
vector dependent variable /. This should come as no 
surprise since the generality of the surrogate- 
equation technique should not prevent the specific 
nature of the particular partial differential equations 
to which it is applied from manifesting itself. Indeed, 
in this case the appearance of the Laplace equation is 
a consequence of the fact that holomorphic functions 
are harmonic. 


Vx(^>y) = Q\^) 

<Pxiix>y) = ^2(y) 

<Pyix,0) = q3(x) 

‘Pyix.ly) = qA(x) 

and also to the additional constraint that the line 

integral of (n = 1 4) around the boundary of 

D must vanish (ref. 19). The solution to this problem 
is unique to within an arbitrary constant and, given a 
solution for (p, the unknowns u and v can be found 
by differentiating <p with respect to x and y, 
respectively. 

An application of the surrogate-equation 
technique, on the other hand, requires the solution of 
a mixed boundary value problem on D. As we now 
deal with the Laplace equation (7) in the two- 
component vector dependent variable /, we require 
that one condition on /be specified at each point of 
the boundary of domain D. Half of the required 
conditions can be immediately obtained from the 
boundary conditions (8), which were applied to the 
original Cauchy-Riemann formulation. If we let f\ 
and denote the first and second components of /, 
respectively, these conditions can be written as; 


/i( 0 ,>-) = gi(y) 
f\{ix>y) = Q2iy) 
/2(x.0) = ^3(x) 
f2{x.ly) = q4ix) 


Problem Specification 

At this point it is instructive to consider an 
application of the surrogate-equation technique to a 
particular boundary value problem for the Cauchy- 
Riemann system of equations (1) and (2). We choose 
the closed rectangular domain D in the x—y plane 
such that 

D=[{x,y) \0<x<lx, 0<y<ly] 

We require that the Cauchy-Riemann system be 
satisfied on D, subject to the boundary conditions 


We are left with the task of formulating the 
remaining four conditions. Since the ultimate 
objective is to obtain a solution to the Cauchy- 
Riemann system on D, it is natural to invoke 
equations (1) and (2), as required, to obtain the 
additional four conditions on the components of /. 
Note that, if we proceed in this manner, no 
additional information that might have an 
overconstraining effect is introduced into the 
problem. In particular, the use of the Cauchy- 
Riemann conditions merely restates the fact that u 
and V are harmonic conjugates. Hence, in the case at 
hand, the remaining conditions are found to be: 


u{0,y) = qx{y) 

u{lx,y) = q2^) 

v{x,0) = q^ix) 
v{x,ly) = qi{x) 


( 8 ) 


Following the classical approach, this Dirichlet 
problem for the Cauchy-Riemann system would be 
transformed into a Neumann problem for the 
potential equation ^Pxx'^'Pyy—^y subject to the 
boundary conditions 






dy 


<7l 


(0,3') 

a 

dx^^ 


dy 


A 


(x,0) 


(x,/v) 
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If we now solve the resultant mixed boundary 
value problem, we will find the two conjugate 
harmonic functions u and v that satisfy the boundary 
conditions (8) imposed on the original Cauchy- 
Riemann system. Hence, we will have found a 
solution to the original Cauchy-Riemann problem by 
solving a second-order partial differential equation 
and without having made use of either the 
irrotationality or incompressibility conditions as field 
equations in our development. This is precisely the 
objective of the surrogate-equation technique: to 
find a second-order partial differential system that 
can serve during the solution procedure as a 
surrogate for the original first-order system while 
neither broadening nor restricting its set of 
admissible solutions. 

For convenience the three boundary value 
problems discussed thus far are schematically 
illustrated in figure 1. 

Computational Results 

Given the potential and surrogate-equation 
formulations of the Cauchy-Riemann problem, as 
described in preceding sections, it is a quite 
straightforward matter to compute approximate 
numerical solutions to both of these boundary value 
problems. 

Such illustrative computations are reported in 
detail in reference 20. There each problem is 
discretized by using second-order accurate central 
differencing, and the resulting systems of algebraic 
equations are solved by successive line overrelaxation 
(SLOR). The problem specification is completed by 


choosing the functions q\ through 94 for Q<x<lx 
and 0<y<ly such that 


9i0) = 


0.5y 


92(y)= 1.0 


, , 0.5x(1.0-0.5A://;f) 

— 

‘y 

, , 0.5xt(1.0-0.5a://v) 

Q4(x)= ^ ^ 

‘y 


0.75 ly 
lx 


The exact solution to this problem is 


0.5>>(1.0-X//;f) X 
u(x,y) = ^ — 


Iv 


ly 


v(x.y) = 


'y ‘X 

0.5x(1.0-0.5x/lx) 
ly 


y(l.0-0.25y/ly) 

lx 

For both the computations with the potential 
formulation and those with the surrogate-equation 
formulation, the SLOR procedure was continued 
until the resultant approximations to the u and v 
velocity components differed from the exact solution 
by at most 1.0x10“® at any point in the 



u ' qj (y) 

/ 


V = Q4 (x) 



V ■ (x) 

Cauchy-Riemann formulation 


u = q2(y) 

\ 


(Py = Q4 (X) 


d , 6 

^ ^ '’4 

?2 = ^4 (X) 


<Px = qj (y) 

tPxx 

= '12 ’y* 

<1 ' Pi fy) 

d 

^xx ’’’ ^yy " ® 




dj( h ‘ dy 'll 



¥>y = Q3 (X) 

Potential formulation 


fl = q2 (y) 

d , d 

'2 =-§-^2 


A, 

dy 'r <13 
f2 = q3 (x) 

Surrogate-equation formulation 


Figure 1. - Three equivalent formulations of the same boundary value problem. 
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computational domain. With the parameters lx and 
ly set at values of 2.0 and 1.0, respectively, 
computations were carried out on three meshes of 
successively doubling point density. 

As expected, the numerical solution to the 
potential formulation behaved such that the u and v 
velocity components, as computed from the x and y 
derivatives of <p, respectively, converged to within the 
above specified tolerance of the exact solution in an 
orderly fashion. Hence the principal nontrivial result 
is that, as expected from the foregoing theoretical 
considerations, the numerical solution to the 
surrogate-equation formulation also behaved well, 
producing the desired approximate solution to the 
required degree of accuracy. 

Summation 

By using the Cauchy-Riemann equations 
describing the two-dimensional, incompressible, 
irrotational flow to provide a model problem, the 
essential features of the surrogate-equation technique 
have been illustrated. The objective of constructing a 
second-order partial differential system to serve 
during the solution procedure as a surrogate for the 
original first-order system while neither broadening 
nor restricting its set of admissible solutions appears 
on the basis of both theoretical considerations and 
the results of computational experimentation to have 
been achieved. 


Euler Problem 

Having illustrated the essential features of the 
surrogate-equation technique by means of the simple 
Cauchy-Riemann model problem, we now proceed to 
the consideration of more realistic and more 
complicated problems, where the utility and 
generality of the surrogate-equation technique can be 
more thoroughly displayed. This section treats the 
computational simulation of the steady flow of a 
perfect fluid under either purely subsonic or purely 
supersonic conditions. Transonic flows are discussed 
in the following section. 

The conventional approach to the computational 
solution of steady, subsonic inviscid flow problems is 
to make use of either the potential or the stream 
function. By doing so, the Euler system of partial 
differential equations, which has only first-order 
derivatives of primitive variables, is replaced with a 
second-order partial differential equation in a 
derived dependent variable. Given this second-order 
equation, one can then draw on the large body of 
experience concerning the design of relaxation 
procedures for such equations in order to arrive at a 


solution algorithm. This convenience is compensated 
for by a loss of generality. The potentitil function 
formulation is limited in application, by definition, 
to irrotational flows. The stream function 
formulation is essentially two-dimensional and 
furthermore for transonic flows is hampered by 
density being a double-valued function of the mass 
flow parameter and by the saddle point that exists at 
the sonic line (refs. 1, 21, and 22). 

An alternative to this approach is to seek a steady 
solution that is the temporal asymptote of solutions 
to the unsteady Euler equations. Assuming that such 
an asymptote exists, this method has the advantages 
both of avoiding the restrictions inherent in the 
potential and stream function formulations and of 
allowing one to deal with purely hyperbolic first- 
order partial differential equations. For such 
equations one may once again draw on a large body 
of experience when designing a solution procedure. 

In the case of time-accurate, time-asymptotie 
solution of the unsteady Euler equations the 
principal disadvantage lies in the long computing 
times to be expected. As a result of efficiency 
comparisons of time-accurate methods with 
relaxation methods for certain model problems, 
Lomax and Steger (ref. 2) found “. . . that relaxation 
methods converge from one to two orders of 
magnitude faster than time-accurate ones.” This led 
them to conclude that “. . . there is a real need to 
improve our relaxation techniques for many types of 
equations modeling steady-state fluid flows.” More 
specifically, they found that, although . . some 
techniques for relaxing isentropic-flow equations in 
terms of the primitive variables have been 
developed” (refs. 1 and 23 to 25), “. . . a suitable 
relaxation procedure for the general Eulerian 
equations has not emerged, so far as the authors 
know” and further that “a major problem area 
where relaxation schemes have yet to be exploited is 
in the calculation of inviscid rotational and energy- 
input flows.” 

The absence of suitable relaxation procedures for 
the steady Euler equations can be accounted for in 
the observation that for regions of subsonic flow, 
where the equations exhibit elliptic behavior, the 
natural differencing techniques, when applied to 
these first-order equations, lead to unstable finite 
difference equations. At this writing, several 
attempts by various researchers at resolving this 
difficulty have, as discussed previously, provided 
interesting results and useful insights but have met 
with less than complete success. 

The realization that the impasse in the 
development of relaxation methods for the steady 
Euler equations is caused by their first-order nature 
leads one quite naturally to consider an application 
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of the surrogate-equation technique. Since this 
technique provides a steady, second-order system of 
partial differential equations whose solution also 
satisfies the Euler equations, we can at once avoid the 
restrictions of the second-order potential and stream 
function formulations and use natural differencing 
techniques to create stable finite difference 
equations. Having thus devised a method for 
obtaining stable, discrete algebraic analogs to the 
steady Euler partial differential equations, we can 
then proceed to investigate various relaxation 
procedures with the aim of identifying those that 
provide a converged solution more efficiently than 
does the time-accurate solution of the unsteady 
equations. Note here that it is not the intention of the 
present work to perform such an investigation of the 
various possible relaxation procedures. Rather, these 
computations are of an illustrative nature. The goal 
then is to present the surrogate-equation technique 
and to show how its application to the equations of 
inviscid fluid flow leads to more powerful and 
general computational procedures than are presently 
available. In particular, we will presently illustrate its 
use for the computation of rotational subsonic flows 
by using the full Euler equations. In the next section 
the surrogate-equation technique is used to create 
stable, fully conservative, type-dependent finite 
difference equations for the numerical solution of an 
inviscid transonic flow problem. 

The purely supersonic flow problem can be readily 
solved by any number of proven techniques. 
Consequently it is included here only to provide a test 
case for the surrogate-equation technique. 


g = 


pv 

puv 

pv^ +p 

(E+p)v 


(10b) 


Here the density, static pressure, and velocity 
components in the x and y directions and the total 
energy per unit volume are denoted by p, p, u, v and 
E, respectively. Furthermore the total energy per unit 
volume can be expressed as 

where the specific internal energy e is related to the 
pressure and density by the simple gas law 


p = (y-l)pe 


with 7 denoting the ratio of specific heats. 

For definitude we now assume that we wish to treat 
a flow that is approximately aligned with the x 
direction and hence wish to rewrite equations (9) 
solely in terms of /. To this end, note (ref. 26) that / 
and g are homogeneous functions of first degree in 
the components of the vector of conservative 
variables w, where w is defined such that 


Equations of Motion 

As is well known, the two-dimensional flow of a 
perfect fluid can be described by specifying four 
partial differential equations, together with the 
appropriate auxiliary relations and boundary 
conditions. These partial differential equations are 
known as the Euler equations and are statements 
concerning the conservation of fluid mass, 
momentum, and energy. For steady flow they can be 
written in vector form as 


fx+gy - 


( 9 ) 


where x and y are the coordinates of a Cartesian 
reference frame. 


/= 



+P 


puv 

{E+p)u 


(10a) 


H' = (p, pu, pv, 


Hence it follows from Euler’s theorem on 
homogeneous functions (ref. 27) that f=Aw and 
g = Bw, where the Jacobian matrices A and B are 
defined such that A =df!dw and B = dg/dw and their 
elements are given explicitly in appendix A. Thus, 
wherever A ~ ’ exists, we can write 


g^Tf 


where T = BA ~ * . Since A is singular only when the 
absolute value of the u velocity component either 
vanishes or is equal to the local sonic velocity 
(appendix A), assuming A to be nonsingular 
introduces no essential limitations for the purely 
subsonic and supersonic flow cases presently under 
consideration. We can now rewrite equation (9) as 


and 


/x + (7/)^=0 
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This equation can in turn be simplified (appendix B) 
to yield 

fx-^Tfy=0 ( 11 ) 

Thus far, we have done no more than to rewrite the 
Euler equations (9) in the form of equation (11), so 
that they are expressed in terms of the single vector of 
dependent variables /. Equation (11), together with 
the flow tangency condition and the appropriate 
upstream and downstream (or far field, if relevant) 
conditions, constitutes a completely specified partial 
differential problem. 

Second-Order System 

The Euler equations, as given in equation (1 1), can 
be rewritten by using differential operator notation 
to yield 


contain those solutions to the Euler equations that we 
seek. Any possible additional elements of the 
solution set to equation (14) that do not satisfy the 
Euler equations are eliminated by means of the 
boundary condition specification, which is discussed 
in the next subsection. 

This surrogate second-order partial differential 
system possesses some interesting properties. By 
virtue of the form of the operator (13) the second- 
order system has no cross-derivative terms. This 
convenience results in a considerable simplification 
of its finite-difference equation analog. Furthermore 
the choice of the operator leads to a pleasing 
behavior of the characteristic directions associated 
with equation (14). We first recall that the 
characteristic directions of the Euler equations are 
determined by the eigenvalues of the matrix A~^B. 
We then observe that the characteristic directions 
associated with the surrogate second-order system 
are determined by the square roots of the eigenvalues 
of the matrix 7^. However, the matrices A~^B and 
T are similar since 


Then, in a manner similar to that used previously for 
the Cauchy-Riemann equations, we can create the 
desired second-order system by operating on 
equation (12) with the differential operator 


d 

dx 



(13) 


to yield 


dx^ 





f=0 


This equation can be reexpressed (appendix B) as 


which can in turn be simplified to 



(14) 


Equation (14) is the surrogate second-order equation 
that will be used here to obtain solutions to the Euler 
equations. Since this surrogate equation has been 
obtained from the original equation by a process of 
differentiation, we expect that its solution set will 


A-^B^A~HT)A 


and hence have the same eigenvalues. It then follows 
that the characteristic directions of the system of 
second-order partial differential equations described 
by equation (14) have slopes equal to ±X, where the 
X are the eigenvalues of the matrix A~^B. This 
means that, in applying the surrogate-equation 
technique to the Euler equations to obtain equation 
(14), we have retained the original Euler 
characteristic directions and added to them their 
reflections through the x axis. For subsonic flows this 
property may not be of great importance. However, 
for supersonic flow or for the extension of the 
method presented here to transonic flow the behavior 
of the characteristic directions gains considerably in 
significance. It is then quite useful to observe that, if 
the coordinate system is chosen such that the flow is 
aligned with the x direction, the standard successive- 
line-overrelaxation iterative procedure possesses the 
same symmetry with respect to the x axis as do the 
characteristic directions of the surrogate second- 
order system. From this point of view, one would 
then expect the introduction of the additional 
characteristic directions associated with the second- 
order system to have a minimal effect on the 
behavior of the SLOR solution procedure. 

Boundary Conditions 

Here we are concerned, in general, with two types 
of boundaries to the physical domain of interest. The 
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first type, which is referred to as a solid boundary, 
occurs at the interface of the fluid with some 
substantial obstacle, such as the wall of a passage or 
the surface of an airfoil. The second type, which is 
called a flow boundary, occurs when for practical 
reasons one arbitrarily prescribes a boundary in the 
fluid-flow field itself that is not a solid boundary but 
beyond which the flow simulation will not proceed. 

We seek a solution to the Euler equations and, as is 
well known (ref. 28), it is a both necessary and 
sufficient solid boundary condition for these 
equations that the fluid-flow velocity vector be 
parallel to the wall slope at the point of contact. 
Hence we require that the surrogate second-order 
system also satisfy this wall tangency condition. 
Furthermore, to insure that we admit no solutions to 
the second-order system that do not satisfy the Euler 
equations, we require that precisely these Euler 
equations also serve as boundary conditions. 
Although no mathematical proof of the sufficiency 
of these additional conditions is presented herein, the 
computational results to be reported serve as a strong 
indication that this is in fact the case. 

Since the flows to be treated subsequently are 
internal flows, the flow boundaries are of the inflow 
or outflow type rather than the far-field type 
associated with external flows. These far-field 
boundaries are normally treated by assuming that 
either “free stream” conditions can be used for the 
values of the unknown or that some far-field 
solution, obtained by other means, is available to 
determine their values. 

The appropriate inflow and outflow boundary 
conditions for internal flows are dependent on the 
exact physical and mathematical nature of the 
problem to be solved. Quite often the inflow 
boundary is treated by specifying the values of all 
unknowns along its extent. Although such 
specification precludes any influence of the flow 
conditions downstream on those at the entrance, this 
nevertheless often results in a physically realistic 
boundary condition. For either supersonic or 
subsonic outflow a commonly used boundary 
treatment, again assuming that it is compatible with 
the physics and mathematics of the flow, is to assume 
that the values of the unknowns do not vary in the 
flow direction. Use is made of such flow boundary 
conditions in the calculations to be discussed 
subsequently. However, the simplified inflow and 
outflow conditions applied successfully to these 
model problems cannot in general be expected to 
yield good results in more complicated flow 
simulations. 

Problem Specification 

To test the surrogate-equation formulation 


described earlier, we consider the simulation of a 
number of two-dimensional internal flows. Since, as 
mentioned previously, very little appears to be 
known about the design of relaxation procedures for 
the Euler equations, the rational development of 
such a procedure requires that simple tests be made 
of the validity of the concepts involved. 

We first compute the flow in the supersonic region 
of a two-dimensional hyperbolic nozzle. Any 
problems due to the introduction of additional 
characteristic directions by the second-order system 
should be revealed by such a case. Furthermore the 
necessary inflow boundary conditions and a series 
solution that is valid close to the sonic line are 
available from the work of Hall (ref. 29). A more 
complete specification of this test case is given in 
reference 30. 

The next test is to compute the purely subsonic 
flow in a two-dimensional symmetric nozzle with 
sinusoidally shaped walls. This geometry is of 
interest because the subcritical flow through it should 
have two axes of symmetry: the nozzle centerline and 
the geometric throat. Further details concerning this 
case are given in reference 3 1 . 

We then consider the computational simulation of 
inviscid shear flows through curved passages. To this 
end, use is made of a class of inviscid, 
incompressible, rotational flows, presented by 
Shercliff (ref. 32), that can be described by the 
stream function 


4/ = C exp (ky) cosh (lx) 


where C, k, and / are constants and x and y are 
Cartesian coordinates. 

As is explained in greater detail in reference 32, 
this stream function describes the flow of 
an incompressible fluid through a bend of angle 
2 arc tan (l/k) that transitions between two 
asymptotic flows that are rectilinear shear flows. The 
bend is quite abrupt and the streamlines approach 
their asymptotes with exponential vigor. Also,as the 
streamtube cross section is greatest at the symmetry 
axis of the bend, these flows first decelerate and then 
reaccelerate as they complete their passage through 
the bend. 

For precision, we have chosen here to examine 
computationally the compressible inviscid flow 
through finite sections of 90° Shercliff bends. The 
results can then be compared with one another and 
for reference with what is subsequently referred to as 
the augmented incompressible flow solution. 

This augmented incompressible flow solution 
consists simply of velocity components obtained 
directly from the stream function that represents the 
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corresponding incompressible flow, together with 
density and static pressure values. These values are 
estimated by specifying density and static pressure 
profiles at the bend entrance and then obtaining their 
distributions throughout the flow field by making use 
of the constancy of both entropy and total enthalpy 
along streamlines while assuming that the perfect -gas 
law applies. In this way we obtain information 
suitable for specifying the entrance conditions for the 
compressible flow computations, for generating 
starting conditions for the iterative solution 
procedure, and for use as a baseline against which we 
can compare the results of our surrogate-equation 
compressible flow computations. 

For variety a number of computations have been 
performed in a different sort of bend, referred to 
herein as a circular-arc bend. This is a bend that 
transitions between two rectilinear sections by 
utilizing a section whose walls are arcs of two 
concentric circles, with tangency being required at 


the joints. Further specifications for these cases are 
given in reference 20. For convenience, all four test 
cases are schematically illustrated in figure 2. 

Computational Results 

Extensive test computations have been carried out 
on the surrogate-equation representation of the Euler 
equations described previously. Of particular interest 
are the results obtained in the four test cases shown in 
figure 2. In all computations reported herein, a 
sheared coordinate system was chosen for the 
physical domain. Such a coordinate system although 
nonorthogonal is simple and convenient and, since 
the bounding walls are coordinate lines, it facilitates 
the precise application of wall boundary conditions. 

In each of the subsonic test cases the governing 
partial differential equations are elliptic and were 
discretized by means of second-order accurate central 
differences. The resulting algebraic equations were 
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(a) Supersonic nozzle (see refs. 29 and 30 for details). (bl Subsonic nozzle (see ref. 31 for details). 




(c) Shercliff bend flow, decelerating flow with turn- (d) Circular-arc bend flow, constant-cross-section 
ing angles up to 9CP (see refs. 20 and 32 for details). flow with turning angles up to 9(f (see ref. 20 

for details). 

Figure 2. - Nomenclature for Euler equation test cases. 
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then solved as a coupled system by using an SLOR 
iteration scheme and a block-tridiagonal analog to 
the Thomas algorithm (ref. 33). Further details are 
contained in references 20 and 30. 

The governing equations for the supersonic test 
case are of course hyperbolic. They were discretized 
by using three-point upwind differences in the flow 
direction and three-point central differences in the 
transverse direction. The resulting coupled system of 
algebraic equations was then solved by using an 
implicit marching scheme. More information 
concerning the details of both the discretization and 
the solution procedure are to be found in refer- 
ence 30. 

For all results presented in this section the 
following normalizations have been employed: 

(1) For velocity components, c., the critical speed 
at the location of maximum entrance velocity 

(2) For density, po, the density at this same 
location 

(3) For static pressure, pqc? 

(4) For length, bend cross section at the axis of 
symmetry or the nozzle throat cross section 

The results of the supersonic nozzle test are 
summarized in figures 3 through 6. There the velocity 
components, density, and static pressure calculated 
by using the surrogate-equation technique are 
compared with the predictions of the Hall solution. 
The agreement is excellent close to the sonic line, 
where the Hall solution is valid. As might be 
expected, the present solution deviates from the Hall 



Figures. - Nondimensional u-velocity distribution - 
supersonic flow. 



Figure 4. - Nondimensional v-velocity distribution - 
supersonic flow. 



Figure 5. - Nondimensional density distribution - 
supersonic flow. 
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subsonic flow. Figure 11. - Nondimensional u-velocity profiles - subsonic flow. 
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solution as the location becomes more remote from 
the sonic line and the nozzle axis of symmetry. No 
detrimental effects attributable to the additional 
characteristic directions introduced by the surrogate- 
equation formulation are apparent in the 
computational results. 

The subsonic nozzle test also produced 
encouraging results. Figures 7 through 10 show the 
computed distributions of velocity components, 
density, tmd static pressure; figure 1 1 presents the u 
velocity profiles obtained at several transverse 
locations. 

To illustrate the utility of the surrogate-equation 
technique for obtaining solutions to the Euler 
equations in flow situations of contemporary 
interest, a computational study has been made of the 
effects of compressibility on the flow through a 
Shercliff 90° bend. The section of the bend used in 
the study, together with the computational grid 
employed, is illustrated in figure 12. The density and 
the static pressure were assumed to hold their normal 
atmospheric values at the bend entrance. The ratio of 
maximum to minimum entrance velocities was set 
equal to 1.5 for the entire study, while the minimum 
entrance velocity was increased from an initial value 
of 50 meters per second through 100, 150, 200, and 
250 meters per second to a final value of 275 meters 
per second. 

While the complete results are given in reference 
20, a sampling of the results of the study is shown in 
figures 12 and 13. There, for the 275-meter -per- 
second case the augmented incompressible flow 
solution and the compressible flow solution obtained 
by means of the surrogate-equation technique are 
shown. The intervening cases, which are not 
illustrated here, show a gradual transition from 
essentially incompressible behavior to the highly 
compressible case of figure 13. 

Results for a 30° circular -arc bend with a 
minimum entrance velocity of 50 meters per second, 
a ratio of maximum to minimum entrance velocities 
of 1.0, and atmospheric values of density and 
pressure at the upstream domain boundary are 
shown in figures 14 and 15. Figure 14 shows the 
solution obtained on the basic computational mesh; 
figure 15 shows the corresponding solution obtained 
on the refined mesh illustrated there. Only the 
relatively minor adjustments in the solution to be 
expected as a result of such a mesh refinement are in 
fact observed. It is, however, interesting to see that in 
both computations a slightly anomalous behavior is 
apparent in the density field. A more refined 
treatment of the wall boundary conditions must be 
considered as a strong candidate for the eventual 
elimination of this behavior. 

During the course of the computations reported 
herein, some data regarding the efficiency, accuracy, 


and stability of the surrogate-equation algorithm for 
the Euler equations have also been collected. 

The algorithm requires approximately 0.0026 
central arithmetic unit (CAU) seconds per grid point 
per iteration when executed using double-precision 
arithmetic on the NASA Lewis Univac 1100/40 
computer system. 

For a simple model problem an experimental 
determination of the actual order of accuracy of the 
surrogate-equation algorithm was performed. For 
fixed flow conditions and a computational domain of 
fixed dimensions, computations were performed on 
three grids, with successive grids having the grid 
point spacings in each direction halved from their 
previous values to yield normalized mesh sizes of 1.0, 
0.5 and 0.25, respectively. The results of this study 
are presented in figure 16, where the error, 
representing the difference between the exact and 
approximate computed solutions at a typical grid 
point, is plotted as a function of the normalized mesh 
size. By virtue of the logarithmic scales used in figure 
16, the observed order of accuracy for the surrogate- 
equation algorithm used in this section can be easily 
estimated to be approximately 2. 

As of this writing the observations made 
concerning the stability of the surrogate-equation 
algorithm are rather qualitative in nature. Formal 
stability bounds have yet to be determined and our 
computational experimentation, although quite 
extensive, does not suffice to estimate them. It does 
appear, however, that by using the surrogate- 
equation technique, we have, as a minimum, left the 
realm of unconditionally unstable finite difference 
analogs to the steady Euler equations and entered 
into the realm of conditionally stable, and apparently 
quite robust, ones. The conditional stability seems to 
be a consequence of the use of the first-order, steady 
Euler differential equations in the formulations of 
the finite difference equations used at solid 
boundaries. Given the rather perverse behavior of 
finite difference analogs to these equations in regions 
where they exhibit elliptic nature, it is quite plausible 
that their introduction would effect such a stability 
reduction. Hence an obvious area for further study is 
the more exact determination of, and subsequent 
improvement in, the stability properties of our 
algorithm. 

Conservation Form 

For the computations discussed in this section, we 
have applied the surrogate-equation technique to a 
nonconservation form, equation (11), of the Euler 
equations. Since the use of this nonconservation 
form presents no serious difficulties for the sort of 
computations discussed here, it was adopted because 
of its simplicity. Should one, however, wish to 
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Figure 14. - Circular-arc bend flow - converged solution, basic mesh. 

differenced by employing the usual techniques (i.e., computational simulations of the steady subsonic or 

second-order, accurate central differencing for supersonic flow of a perfect fluid. Solutions to the 

purely subcritical flows) to yield fully conservative Euler equations are obtained without resort to either 

(refs. 13 and 34) finite difference equations. the potential or stream function formulations, and 

the inherent limitations of these formulations are not 
Summation shared by the present approach. Furthermore the 

time-asymptotic solution of the unsteady Euler 
In this section we have shown that the surrogate- equations is also avoided, and convergence 

equation technique can be used to perform acceleration by either relaxation or some other non- 
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time-accurate process is possible. 

The algorithm presented here is reasonably 
efficient but possesses ample possibilities for 
improvement. It is approximately second-order 
accurate and is conditionally stable, with further 
study being necessary to determine precise stability 


bounds. Although the results presented here were 
based on a nonconservation form of the Euler 
equations, the development of an algorithm that uses 
a conservation form, as discussed in this section, 
presents no essential difficulties. 
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Figure 16. - Experimental accuracy determination for Euler 
problem. 


Transonic Flow 


reach of these methods. 

There is a need for a technique for computa- 
tionally simulating transonic flows that is, on the one 
hand, computationally efficient and, on the other 
hand, not subject to the restrictions of a potential 
flow formulation. In this section, we demonstrate the 
ability of the surrogate-equation technique to serve 
this need. 

For simplicity and clarity the surrogate-equation 
technique is first applied within the context of 
transonic small-disturbance theory, where 
comparison can be made with the conventional 
perturbation-potential approach. In this manner the 
validation of the method for use in transonic flow is 
separated from the demonstration of its usefulness in 
rotational flow, which was presented in the preceding 
section. Its extension to the full steady Euler 
equations is discussed later in this section. 


-Applications of the surrogate-equation technique 
are by no means limited to the classes of steady 
subsonic and supersonic flows from which the 
examples discussed in the previous section were 
drawn. In fact, as illustrated here, the technique can 
be used to good advantage for the computational 
simulation of steady transonic flows. Such flows are 
of considerable mathematical interest and have great 
technological importance. 

As discussed previously the main body of methods 
for the computation of steady inviscid transonic flow 
that have been developed thus far either utilize the 
time-asymptotic solution of the unsteady Euler 
equations to obtain a steady solution or iteratively 
sohe the steady potential or perturbation potential 
equation. In the former case computation times are 
long; in the latter case the generality of the model 
equations is sacrificed in the name of computational 
efficiency. 

The transonic small-disturbance theory, upon 
which the perturbation-potential approach is based, 
is derived under the assumptions that body surface 
slopes are everywhere small (so that flow quantities 
are small perturbations about their free-stream 
values) and that the free-stream Mach number is near 
unity. In practical situations these assumptions are 
not always strictly met, but nonetheless many cases 
of engineering interest can be adequately treated. 
Where the assumptions of small-disturbance theory 
are grossly violated, resort is made to the full 
potential formulation. 

The assumption of potential flow inherent in both 
of these approaches and the resultant restrictions to 
irrotational and isentropic conditions can prove quite 
troublesome under certain circumstances. Clearly 
shear flows cannot be treated, and hence a large class 
of technologically important flows lies beyond the 


Small-Disturbance Approximation 

The details of small-disturbance theory have been 
thoroughly discussed in the literature (refs. 35 to 39), 
and our remarks here are accordingly limited. In 
brief, if we assume that the flow of interest can be 
represented as a disturbance to a uniform flow, that 
in particular the disturbance velocity components are 
small with respect to the mean velocity, and that for 
transonic flow the Mach number of the undisturbed 
flow is close to unity, the exact Euler equations can 
be replaced by approximate transonic small- 
disturbance equations. These equations, although 
quite simple, retain the essential nonlinear, mixed 
elliptic-hyperbolic character of the exact equations. 
Furthermore weak solutions of the transonic small- 
disturbance equations that admit discontinuous 
jumps approximating shock waves can be obtained. 

Although alternative forms of this equation are 
readily available and have in fact been employed 
successfully by other investigators, the following 
formulation is sufficient for our present purposes: 

l-Mi-(7+l)Mi«j^ + g=0 (17) 

Flere u and v represent the disturbance velocity 
components in the x and y directions, respectively, 
normalized by the uniform stream velocity, referred 
to as Ua,. Also, Moo is the Mach number of the 
undisturbed uniform stream and y is the ratio of 
specific heats. Equation (17) can be rewritten in 
conservation form as 


ax ay 


(18) 
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T 


where 


We can supplement this single equation in two 
unknowns by the irrotationality condition 


dv du 
dx By 


(19) 


which is also a consequence of the small-disturbance 
assumption, to obtain a closed set of equations. 
Equations (18) and (19) then constitute a system of 
first-order, partial differential equations of mixed 
elliptic-hyperbolic type. They are elliptic or 
hyperbolic according to whether the u component of 
the disturbance velocity is smaller than or larger than 
the critical disturbance velocity «*, which is defined 
as 


(T+l)Mi 


Perturbation-Potential Equation 

Although a computational procedure could be 
devised to calculate an approximate solution to the 
system of equations (18) and (19) directly, this is not 
normally attempted. This is so because equation (19) 
implies the existence of a scalar perturbation velocity 
potential >f> such that ^x = u and tpy = v. Hence 
equation (18) can be rewritten as 

|^(l - Mi - ^ Mi ^ J = 0 (20) 

while equation (19) reduces to the identity 
'Pyx ~ ‘Pxy ~ 0 

Equation (20) is a second-order partial differential 
equation of mixed elliptic-hyperbolic type in the 
scalar unknown <p. This appears to be advantageous 
on two counts. First, in solving equation (20) for >p 
and thence for u and v we can possibly obtain these 
velocity components with less computational effort 
than would be required by a direct solution of the 
system of equations (18) and (19). Second, since 
equation (20) is of second order, the formulation of 
stable, conservative finite difference equations for its 
discrete representation is greatly facilitated. 


Although the transonic small-disturbance 
equations have been known for a considerable time, 
their essential nonlinearity has impeded progress on 
their analytical solution. On the other hand, their 
mixed elliptic-hyperbolic nature, together with the 
fact that the locations at which type changes occur 
cannot generally be prescribed apriori, confounded 
attempts at their approximate numerical solution. In 
1969 Cole (ref. 35) summarized the status of 
transonic small-disturbance theory and set up the 
problem of plane mixed flow past an airfoil, 
including a discussion of the far field. Subsequently 
Murman and Cole (ref. 3) devised a numerical 
method for the computation of an approximate 
solution to this problem. The details of their basic 
method and its subsequent improvements and 
generalizations are given in the literature, 
particularly in references 3 and 40 through 44. A 
brief discussion of those aspects of the method that 
are germane to our present purpose follows. 

Murman and Cole overcame the difficulties 
associated with the mixed elliptic-hyperbolic nature 
of the transonic perturbation -potential equation by 
introducing the idea of type-dependent differencing. 
By automatically adapting the finite difference 
equations at each grid point of the computational 
domain to suit the local nature of the flow, they were 
able to construct an iterative procedure for the 
solution of mixed flow problems that “captures” any 
shocks that may be present and represents them as 
steep gradients and that is computationally stable 
and in conservation form. 

For the field equations and boundary conditions, 
second-order accurate central differencing is used for 
derivatives in the y direction and for derivatives in 
the X direction in regions where the flow is subsonic. 
Backward (or upwind) differencing of either first- or 
second-order accuracy is used in the x direction in 
regions where the flow is supersonic. At “sonic” and 
“shock” points special switching operators that 
preserve the conservative nature of the differencing 
scheme are employed. Hence domains of dependence 
are everywhere respected and intercellular fluxes are 
properly conserved. A further consequence of the 
conservation form of the partial differential and 
finite difference equations is that, should shocks be 
present, the proper (isentropic) jump conditions are 
attained. 

The finite difference equations are written in 
implicit form, thereby avoiding the restriction of 
vanishingly small mesh width in the x direction upon 
approaching sonic velocity, which would be 
encountered with an explicit formulation. This 
system of algebraic equations is then solved 
iteratively by the method of successive line 
relaxation. In this fashion the approximate numerical 
solution is recomputed along lines transverse to the 
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flow direction as the computational domain is 
repeatedly traversed in the direction of the flow. In 
subsonic regions the solution is overrelaxed to 
accelerate convergence. 

The boundary condition treatment is such that 
Neumann conditions on the perturbation potential 
are specified on the surface of an immersed body and 
applied, as is consistent with the small-disturbance 
approximation, along a coordinate line. At some 
finite distance from the body a domain outer 
boundary is chosen along which a far-field solution is 
used to provide the necessary boundary condition. 

As singularities can be present at the leading or 
trailing edges of an airfoil about which the flow field 
is to be computed, Murman and Cole chose to locate 
the boundary points of their computational domain 
one half of a cell width from the leading and trailing 
edges. Although this might seem to be a rather 
simplistic remedy, the nature of the singularities in 
question is such that this approach is reasonably 
good. Further details concerning the inner boundary 
condition specification are given in the original 
Murman and Cole article (ref. 3) and in particular in 
the work of Krupp (ref. 41). 


Since /3 is not a homogeneous function of first degree 
in M, equation (21) is not equivalent to a conservation 
form of the transonic small -disturbance equations. 
However, as becomes apparent in the subsequent 
discussion, this is not an impediment to the 
formulation of a conservative second-order system. 
In any case this lack of homogeneity is not present in 
the full Euler equations. 

To create a surrogate second-order equation for 
equation (21), we operate on it with a differential 
operator of the form 


to yield 



+ MB^) 




dx ^ dx 

dy^ 

dy 

^ dx 

dy{ 


( 22 ) 


Second-Order System 

Preparatory to our ultimate goal of using the 
surrogate-equation technique to devise an iterative 
scheme for the Euler equations, we illustrate the 
basic process on the simpler, but nevertheless similar, 
transonic small-disturbance equations. In this way, 
we can develop the method, test it, and compare its 
performance with the Murman and Cole approach. 

The transonic small-disturbance equations (18) 
and (19) can be written in vector form as 


For any nonsingular choice of the matrix N, if M is 
chosen such that 

(23) 

the characteristic directions of equation (22) are 
determined by the expression 


9x2 dy^ 



where 




0 1 





1 

0 


(21) which reduces to 



-> 5)2 


ii 

dy^ 


Hence, here again, as was the case in the examples 
considered previously, the set of characteristic 
directions of the surrogate second-order system 
contains those of the original first-order system as a 
subset, with the added characteristics being the 
reflections of the original ones through the x axis. As 
mentioned previously this symmetry in the 
characteristic directions of the second-order system is 
of the same nature as the symmetry inherent in the 
successive-line-relaxation solution procedure to be 
employed here. Therefore one would expect the 
additional characteristic directions associated with 
the second-order system to have a minimal effect on 
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the solution procedure. Furthermore for the choice 
of the matrix M indicated in equation (23) 



so that the surrogate second-order system has no 
cross-derivative terms. This will serve to simplify its 
finite difference representation. The resulting 
second-order system can be written as 



/=0 (24) 


This equation is in conservation form and can be 
differenced to yield fully conservative finite 
difference equations. 

We note in passing that, although the use of 
equation (23) results in several pleasant 
consequences, one may, under different 
circumstances, wish to consider other specifications 
of the matrix M. 

At this point it is instructive to write out the scalar 
equations represented by equation (24) and examine 
them. It is a straightforward matter to perform the 
necessary algebra to obtain 


dx\du dx' dy 2 
and 

dx^du dx' ^ dy2 


(25) 


(26) 


These equations are quite simple and present obvious 
differencing possibilities, as discussed subsequently. 

Having formulated the surrogate second-order 
system for the transonic small -disturbance equations, 
we now proceed with a discussion of the boundary 
conditions necessary to completely specify the partial 
differential problem being considered for numerical 
solution. As we are numerically investigating a two- 
dimensional internal transonic flow, the boundary 
condition discussion is presented in such a context. 

At the upstream and downstream flow boundaries, 
which are located far from any disturbance to the 
flow field and in regions where the velocity is 
uniform, we require that both the u and v 
disturbance velocity components vanish. At solid 
boundaries we require that the v component of the 
disturbance velocity be equal in magnitude to the 


local boundary slope. This is the usual solid-wall 
boundary condition of transonic small-disturbance 
theory. It then remains to specify conditions on the u 
disturbance velocity components at the solid 
boundaries. These are easily obtained from the 
original first-order system. One may, for instance, 
require that at one wall u be such that equation (18) is 
satisfied and at the opposite wall u be such that 
equation (19) is satisfied. 

This then completes the specification of all 
necessary boundary conditions and furthermore does 
so in a manner designed on heuristic grounds, as 
discussed previously, to restrict the admissable 
solutions to our second-order system to be identical 
with those of the original first-order transonic small- 
disturbance equations. 

Having discussed the formulation of the surrogate 
second-order system for the transonic small- 
disturbance equations, we now proceed to define the 
physical problem that will be used as a vehicle for 
testing the efficacy of the surrogate equation 
technique for inviscid transonic flow computation. 


Problem Specification 

Consider an inviscid flow in a two-dimensional 
channel with a uniform inlet velocity (/„ and inlet 
Mach number Moo . The upper surface of a bicircular 
arc airfoil is mounted on the lower channel wall. The 
channel height is one airfoil chord length, and the 
upstream and downstream flow boundaries are 
located one chord length upstream of the airfoil 
leading edge and one chord length downstream of the 
airfoil trailing edge, respectively. The airfoil half- 
thickness is equal to 10 percent of its chord length. 
Alternatively, this problem can be viewed as 
representing the flow past a 20-percent-thick 
bicircular-arc airfoil mounted at zero angle of attack 
on the centerline of a two-dimensional wind tunnel or 
as an unstaggered linear cascade with a gap-to-chord 
ratio of 2. The problem is schematically depicted in 


V = Wall slope 

® + ^ . 0 
dx dy 
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Figure 17, - Transonic flow problem. 
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figure 17, where the lengths lx and ly assume the 
values 3 and 1, respectively. 

We will examine cases where the Mach number 
Moo is low enough so that, although local regions of 
supersonic flow may be present, the passage is not 
choked. Hence, the flow will be globally subsonic 
and in particular subsonic conditions will prevail at 
both the upstream and the downstream boundaries 
of the domain. 


Discrete Formulation 

Second-order systems like equations (25) and (26) 
can be quite easily discretized, in the spirit of 
Murman and Cole, for transonic flow computations. 
One such possible discretization is presented here. 

We seek an approximate solution at a finite 
number of points distributed over the closed 
rectangular computational domain D, where 


^ = [(Jf.3') ko<x<xo+/;c, yo <y<yo+ly] 


We call these points the grid (or mesh) points and 
define them to be the ordered pairs (x/, yj) such that 


= Uo ~ 

yy = (yo- ^y) +j^y 


where /= 1,...,M 
where j=\,...,N 


where bx and 6>' are the mesh spacings in the x and y 
directions, respectively. The value of some dependent 
variable a at the point (x,-, yj) is denoted as a-,j, that 
is. 


Qij =a(Xi,yj) 


Furthermore no notational distinction is made 
between the solution to the partial differential 
equations (25) and (26) and the approximate solution 
to the finite difference equations now to be 
constructed. 

In regions of subsonic flow, where equations (25) 
and (26) exhibit elliptic behavior, second-order 
accurate central differencing can be employed in 
both the X and y directions to yield as finite 
difference representations 




Here the mesh ratio r is defined such that r=bx/by, 
and a represents either u or v. Furthermore, to 
second-order accuracy, we may write 



while d^/du can be expressed as 
|^=l-Mi-(T+l)Mi« 

On the other hand, in regions of supersonic flow, 
where the governing equations exhibit hyperbolic 
behavior, we proceed differently. Although the 
second-order accurate differencing in the y direction 
is retained, in the x direction we switch to backward 
(or upwind) differencing, which is first-order 
accurate. Higher-order accurate differencing could 
be used; however, first-order accuracy may be 
preferable and is in any case sufficient for our 
purposes. This asymmetric differencing is designed 
to respect the domain of dependence of the partial 
differential equations in the hyperbolic region while 
yielding stable finite difference equations. In this 
fashion we obtain 



( 28 ) 


(v 
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Equations (27) and (28) determine the basic nature of 
the discrete approximation to the field equations. 

As is now well-known (refs. 40 and 43), in 
transonic flow computations based on inviscid model 
equations it is not sufficient, in order to insure that 
the proper weak solution be captmed, that the partial 
differential equations be written in conservation 
form and then differenced in a type-dependent 
fashion. One must, at points of transition between 
regions where either one or the other form of the 
difference equations is used, introduce special 
switching approximations. The function of these 
equations is to insure that the intercellular flux 
balances are not disturbed by the change from one 
form of difference operator to another. This in turn 
is intended to assure that the weak solution that is 
captured by the numerical procedure will depend 
only on the original partial differential equations and 
boundary conditions and not on local perturbations 
caused by flux imbalances. 

The form of the switching equations required 
depends on the discretizations used in subsonic and 
supersonic regions and on whether the transition is 
from subsonic to supersonic flow, which is referred 
to as a sonic point, or from supersonic to subsonic 
flow, which is referred to as a shock point. As the 
names imply, sonic points occur at smooth 
transitions, while shock points occur at abrupt ones. 

For the case at hand special attention need only be 
given to the treatment of x derivatives, since no 
change is made in the discretization in the direction. 
Hence for the interior discretization we consider a 
term of the form 

dx \du dx' 

where a represents either the u or the v perturbation 
velocity component. It is a straightforward matter to 
determine that intercellular fluxes will be conserved if 
we discretize this term such that at sonic points 


^\dudx/ .. 

-I i,j 



and at shock points 



For the boundary condition discretization we must 
construct sonic and shock point representations for 
terms of the form da/dx where a represents either v 
or /3- Proceeding as in the previous case, we can 
readily determine that the required representations 
are such that at sonic points 

\dx/.. 2dx 
ij 


and at shock points 


( ^ ) _ ^i+ \,j +g/,,/— 2g/-ij 

\dx^ 28x 

i,j 


To decide in a consistent fashion which of these 
approximations is to be used at a given point {Xj, yj) 
in the computational domain, we employ the 
following switching scheme: 


(^) >0 and 

\ du/ 

(f) and 

\du/ 

1-1.7 

( ^ ) <0 and 

\du/. 

i-Uj 



Then 

Subsonic 

Sonic 

Supersonic 

Shock 


Here dl3/du changes sign from positive to negative as 
u exceeds u*, the critical disturbance velocity. 

Both in the present case and in the case of the 
perturbation-potential equation, the switching 
approximations that are obtained may not be locally 
consistent with the differential equations that they 
model. However, as we seek a weak solution to the 
conservation form of the equations, it is precisely 
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these switching approximations that allow the finite 
difference equations to be consistent in a global or 
integral sense. 

Solution Procedure 

As in most of the previously considered examples 
the finite difference equations are solved by the 
method of successive line relaxation. We sweep 
repeatedly across the domain, moving in the flow 
direction, relaxing the solution on transverse 
(i = constant) lines in order to accelerate the 
convergence of the iterations. The boundary 
conditions are fixed during any given sweep of the 
computational domain and are then recomputed 
after the completion of each sweep. 

The finite difference equations for the totality of 
the interior points along a given transverse mesh line 
can then be written in the form of a matrix equation: 


Mf + ‘ =Ri 


where Af is a block-tridiagonal matrix consisting of 
2x2 matrices, 


af”'*'* = 


A/Tn-^2 


r/z + 1 


•/l + l 


^T 


A/ 


n + \ 
i,N - 1 


Afn + l^fn + l_fn 

■'ij ■'iJ ■' U 


f u — 


Hj 




and Rj is the resulting right-hand-side vector. In 
general a block-tridiagonal matrix elimination 
procedure would be used to solve such an equation. 
However, in the present case the equations for u and 
V are coupled only through their right-hand sides, 
and the matrix elements of the block-tridiagonal 
matrix Mare therefore diagonal matrices. Hence we 
need only solve two simple tridiagonal equations; a 
reduction in the requisite computational effort. We 
accelerate the convergence of the iteration procedure 
such that, if we refer to the r£sult of the tridiagonal 
elimination procedure as the relaxation 

process takes the form 


F," + ‘ =wF7 + ^ -Kl -w)Ff 


To compare the results obtained by using the 
surrogate-equation technique with those of the 
conventional perturbation-potential approach, two 
FORTRAN IV computer programs have been written 
for use on the NASA Lewis Univac 1100/40 
computer system. Both have been written in the same 
spirit, and an equal degree of effort has been 
expended on making each program computationally 
efficient. 

Computational Results 

To generate information that can serve as a basis 
for evaluating the appropriateness and efficacy of the 
surrogate-equation approach to the numerical 
solution of transonic flow problems, two series of 
computations have been performed. For the same 
physical problem, computational domain, initial 
guess, and grid system, computations have been 
carried out for both subcritical and supercritical flow 
conditions with both the surrogate-equation program 
and the perturbation-potential program. Both 
because of the stringent convergence criterion 
employed here and for reasons of general prudence, 
all computations were done using double-precision 
arithmetic. 

Four grid systems were employed during this 
study; extra coarse, coarse, medium, and fine. 
Essentially, each successive grid is constructed by 
halving the grid spacing of its predecessor. More 
detailed information on these grid systems is 
contained in table I. Note that we have sacrificed 
efficiency to simplicity by maintaining a uniform grid 
spacing everywhere rather than stretching the 
streamwise grid spacing upstream and downstream 
of the wall-mounted airfoil. 

We consider one subcritical and one supercritical 
flow case, with upstream Mach numbers of 0.65 and 
0.70, respectively. These two flow cases are described 
in table II. Note that the supercritical case does not 
represent a choked flow. A preliminary survey, using 
both the surrogate-equation program and the 
perturbation-potential program, indicated that 
choking occurs between the upstream Mach numbers 
0.71 and 0.72. This is in reasonable agreement with 
the choking Mach number of 0.73 that is obtained by 
using an approximate relation derived by Spreiter, 
Smith, and Hyett (ref. 45) together with a constant 
for circular-arc profiles that has been estimated by 
Collins and Krupp (ref. 46) on the basis of 
experiments. 

Before proceeding to the more interesting 
supercritical flow results, let us consider the 
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TABLE I.— DEFINITION OF MESHES FOR TRANSONIC FLOW COMPUTATIONS 


Mesh 

Number of 

Number of 

Total 

Mesh spacing 

Mesh spacing 

designation 

points in 

points in 

number of 

in X direction. 

in y direction, 


X direction 

y direction 

points 


6, 

Extra coarse 

30 

6 

180 

0.1 

0.2 

Coarse 

60 

11 

660 

.05 

.1 

Medium 

120 

21 

2520 

.025 

.05 

Fine 

240 

41 

9840 

.0125 

.025 


A 


TABLE II.— NOMENCLATURE FOR TRANSONIC 
FLOW TEST CASES 


Flow 

Upstream 

Critical u 

Critical 

designation 

Mach number 

disturbance 

velocity 

pressure 

coefficient 

Subcritical flow 

0.65 

0.569527 

-1,139053 

Supercritical flow 

.70 

.433673 

- .867347 



Figure 18. - Surface pressure coefficient distribution - 
surrogate-equation formulation, subcritical flow case. 


subcritical flow solutions generated by the two 
programs. Figure 18 presents the pressure coefficient 
distributions on the streamline adjacent to the airfoil 
surface that result from computations on the extra- 
coarse, coarse, medium, and fine grids when using 
the surrogate-equation program. A similar sequence 
of results generated by the perturbation potential 
program is presented in figure 19. In each case the 
solution appears to converge to an asymptote as the 
grid is refined. Furthermore the expected fore-and- 



Figure 19. - Surface pressure coefficient distribution - 
perturbation-potential formulation, subcritical flow case. 


aft symmetry is present in both cases and stagnation 
points are well represented. Figure 20 shows the fine- 
grid pressure coefficient distributions for both the 
perturbation-potential and the surrogate-equation 
formulations. The results agree well with one 
another, an indication that both methods tend to the 
same asymptotic solution. 

We now examine the supercritical flow case where 
a shock is present. In a fashion analogous to that of 
the subcritical flow case figures 21 and 22 present the 
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Figure 20. - Comparison of surface pressure coefficient 
distributions - subcritical flow case, fine grid. 


Figure 22. - Surface pressure coefficient distribution - 
perturbation-potential formulation, supercritical flow case. 



Figure 21. - Surface pressure coefficient distribution - 
surrogate-equation formulation, supercritical flow case. 


pressure coefficient distributions on the streamline 
adjacent to the airfoil surface that result from 
computations on the aforementioned four-grid 
sequence obtained by using the surrogate-equation 
program and the perturbation-potential program, 
respectively. Here again, each solution appears to 
converge to an asymptote as the grid is refined. 

The fine-grid pressure coefficient distributions 
obtained by means of the two approaches are 
compared in figure 23. The results agree well with 
one another. The shock location is reasonably well 



Figure 23. - Comparison of surface pressure coefficient 
distributions - supercritical flow case, fine grid. 


predicted and the shock strengths, although not 
identical, are in fair accord. There are several 
possible explanations for the observed 
underprediction of shock strength and slight 
difference in shock location obtained with the 
surrogate-equation algorithm. For example, the 
specific form chosen for the discretization at shock 
points can influence the behavior of the solution near 
the shock. In particular, although the shock point 
discretization is chosen such that a fully conservative 
difference scheme results, it is not consistent with the 
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governing partial differential equations to the same 
order as the subsonic and supersonic discretizations. 
This consistency problem also manifests itself at 
sonic points, but because of the smooth nature of the 
solution there it causes no extraordinary difficulties. 
Other factors contributing to the discrepancy in 
shock strengths may include differences in the 
boundary condition formulation and implementation 
in the two programs, differing dissipative properties 
due to different truncation errors creating different 
^ artificial viscosities, and the fact that the lack of 

homogeneity of /3 causes equation (21) to be not 
equivalent to a conservation form of the transonic 
small-disturbance equations. The prediction of shock 
strengths by means of the surrogate-equation 
algorithm is presently the object of further study and 
analysis. 

It is interesting that post-shock reexpansion, first 
reported in the experiments of Ackeret, Feldmann, 
and Rott (ref. 47) and later discussed from a 
theoretical point of view by Zierep (ref. 48) and 
Oswatitsch and Zierep (ref. 49), is present in both 
solutions. 

Hence, it appears that, when measured against the 
standard perturbation-potential approach, the 
surrogate-equation technique yields a viable method 
for obtaining numerical solutions to transonic flow 
problems. Recall that the main attraction of the 
surrogate-equation technique rests not on its value in 
providing an alternative method for the solution of 
the transonic small-disturbance equations but rather 
on its applicability to the solution of the full steady 
Euler equations. Nevertheless, a consideration of the 
computational efficiency of the perturbation- 
potential and surrogate-equation programs yields 
some interesting results. 

In the absence of evidence to the contrary, one 
might well expect the perturbation-potential 
formulation to be superior to the surrogate-equation 
formulation when considered from the viewpoint of 
computational efficiency. Such an expectation 


would, no doubt, be based on the observation that 
less work is required to solve the one equation of the 
perturbation-potential formulation than is needed to 
solve the two equations that result from an 
application of the surrogate-equation technique. This 
is in fact the case. We have determined by 
computational experimentation that, in their present 
form, the perturbation-potential and surrogate- 
equation programs require 0.000175 and 0.000256 
CAU second per iteration per grid point, 
respectively. Hence by this measure the surrogate- 
equation program requires approximately 1.46 times 
as much CAU time per iteration per grid point as the 
perturbation-potential program does. 

To determine the computational effort each 
program requires to produce a converged solution, 
we must also take into account the number of 
iterations needed to reduce the residuals to the level 
at which we declare the solution to have been 
reached. To this end, a study involving both the 
subcritical and supercritical flow cases, all four 
computational grids, and both programs has been 
conducted. The results are shown in table III, where 
for each case the number of iterations required under 
optimum relaxation to reduce the absolute value of 
the maximum residual to 10~*<^ is indicated. This 
convergence criterion is rather severe, and a much 
milder one would suffice for engineering 
applications. The use of such a severe convergence 
criterion here can be viewed as a further test of the 
algorithms. 

We note immediately from table III that the 
perturbation-potential program required, for the 
cases shown, on the average 4.36 times as many 
iterations in subcritical flow cases and 5.52 times as 
many iterations in supercritical flow cases as the 
surrogate-equation program. These results are 
illustrated graphically in figures 24 and 25 for the 
subcritical and supercritical flow cases, respectively. 
This marked superiority of the surrogate-equation 
program can be attributed to the fact that most of the 


TABLE in.— CONVERGENCE BEHAVIOR OF TRANSONIC 
FLOW ALGORITHMS 


r 



Number of iterations required for convergence 




Subcritical flow 


Supercritical flow 



Extra- 

coarse 

mesh 

Coarse 

mesh 

Medium 

mesh 

Fine 

mesh 

Extra- 

coarse 

mesh 

Coarse 

mesh 

Medium 

mesh 

Fine 

mesh 


Perlurbalion- 

potential 

formulation 

261 

669 

1879 

6348 

423 

1098 

3774 

12 910 


Surrogate- 

equation 

formulation 

104 

209 

4S7 

838 

130 

312 

639 

1 313 


27 


Formulation 



Figure 24. - Convergence behavior of transonic flow 
algorithms - subcritical flow case. 



Figure 25. - Convergence behavior of transonic fiow aigorithms - 
supercritical flow case. 
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boundary conditions employed in this formulation 
are Dirichlet ones while all of those in the 
perturbation-potential program are Neumann ones. 
It further appears that the convergence of the 
surrogate-equation formulation is less affected by the 
presence of shocks than is that of the perturbation- 
potential formulation. 

In figure 26 the behavior of the absolute value of 


the maximum residual is shown as a function of 
iteration level for both the perturbation -potential 
and surrogate-equation programs. For this 
comparison, we have chosen the supercritical flow 
case on the medium grid. From the data we estimate 
the spectral radii for the iteration matrices of the 
perturbation-potential and surrogate-equation 
programs to be 0.996 and 0.968, respectively. These 
in turn imply asymptotic convergence rates of 
1.74x10“^ and 1.41x10“^. Hence we estimate 
that the surrogate-equation program converges 8.1 
times as fast as the perturbation-potential program. 

By taking account of both the CAU time per 
iteration per grid point and the number of iterations 
necessary for convergence, we estimate that, for the 
cases considered here, the computational effort 
required by the surrogate-equation program varies 
from 0.15 to 0.58 of that required by the 
perturbation-potential program. Thus, even in this 
test case, chosen for simplicity and ease of 
comparison rather than to illustrate the power of the 
surrogate-equation technique, its application to the 
transonic small-disturbance problem results in a 
more efficient algorithm than is obtained by means 
of the conventional perturbation-potential 
formulation. 

As we have no exact solution with which to 
compare the results, we can here only estimate the 
order of accuracy of the programs based on the 
computational results obtained on the sequence of 
four progressively finer grids described in table I. We 
do this by assuming the solution obtained on the fine 
grid to be exact and then calculating the absolute 
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Figure 27. - Order of accuracy estimation - surrogate-equation 
formuiation, subcritical flow case. 
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Figure 28. - Order-of-accuracy estimation - perturbation- 
potential formulation, subcritical flow case. 
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value of the difference between this solution and the 
solution on the various other grids at some 
representative location. This quantity is called the 
representative error and is plotted as a function of 
normalized mesh size for several cases in figures 27 
through 30. 

Comparison of figures 27 and 28 for subcritical 
flow shows that the surrogate-equation algorithm 
appears to be second-order accurate in the 
normalized mesh size but the perturbation-potential 
algorithm approaches first-order accuracy. This can 
be attributed to the fact that in the perturbation- 
potential formulation the flow variables from which 
the measure of error has been calculated must be 
obtained from the computed perturbation potential 
by differentiation. 

Figures 29 and 30 illustrate the results obtained 
from applying this order-of-accuracy estimation to 
the supercritical flow test case. The surrogate- 
equation algorithm retains its second-order accuracy 
in the subsonic region while, here again, the 
perturbation-potential algorithm approaches first- 
order accuracy. We have also attempted to estimate 
the order of accuracy of the two programs at a point 
within the supersonic region. Because of the limited 
size of this region and the large difference in 
resolving power between the coarsest and finest 
meshes, this estimation proved very difficult to carry 
out and its results should be considered to be quite 
approximate. These remarks notwithstanding, the 
results of the order-of-accuracy estimations in the 
supersonic region are also shown in figures 29 and 


30. For each algorithm the order of accuracy in the 
supersonic region is lower than that in the subsonic 
region. This is to be expected since the differencing 
used in the x direction at supersonic points is only of 
first-order accuracy. 

Murman (ref. 43) has stated that based on a linear 
stability analysis the Murman and Cole approach to 
solving the perturbation-potential equation, a 
version of which has been presented here, results in 
unconditionally stable finite difference equations. A 
similar result should hold for the surrogate-equation 
algorithm. Practical confirmation of the 
unconditional stability of the surrogate-equation 
algorithm was provided by the computations 
conducted during the course of the work reported 
herein. 

Application to the Full Euler Equations 

Having established that the surrogate equation 
technique can be successfully applied to the transonic 
small-disturbance equations, we now make use of the 
technique to develop a surrogate second-order system 
for the full, steady, two-dimensional Euler equations 
that will be valid for use in transonic flow 
computations. The resulting second-order system will 
be in conservation form, and fully conservative finite 
difference equations can be constructed for it in a 
fashion analogous to that previously discussed in the 
context of the small -disturbance equations. 

The two-dimensional, steady Euler equations are 
normally written in conservation form as 



Figure 29. - Order-of-accuracy estimation - surrogate-equation 
formuiation, supercritical flow case. 



Figure 30. - Order-of-accuracy estimation - perturbation- 
potential formulation, supercritical flow case. 
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Furthermore as discussed previously /and g are both 
homogeneous functions of first degree in the 
components of the vector w, where w is defined as 
w = [p, pu, pv, E]^. It then follows from Euler’s 
theorem on homogeneous functions (ref. 27) that 
f=Aw and g=Bw, where A and B are the Jacobian 
matrices shown explicitly in appendix A. Hence we 
can reexpress the Euler equations in conservation 
form as 


(Aw)x + {Bw)y =0 

or in operator notation as 


Here again, as in the examples previously discussed, 
the set of characteristic directions of equation (31) 
contains those of the first-order system, equation 
(29), as a subset, with the additional characteristic 
directions being the reflections of the original ones 
through the x axis. 

In the preceding formulation the inverse of the B 
matrix appears. Hence such a formulation would be 
inappropriate at points where the B matrix is 
singular. Such points occur where the v velocity 
component either vanishes or becomes sonic. We can 
avoid this undesirable behavior by making different j 

choices for the matrices A/ and Nin equation (30). If, 
for example, we choose M=A and N=^B, we obtain 




w = 0 


(29) 


We now construct a surrogate second-order system 
by applying the differential operator 


dx 


(M)-|(N) 


to equation (29) to obtain 



d_ 

dx 




Although the cross-derivative terms do not cancel in 
equation (32), it does constitute a usable surrogate 
second-order system for equation (29) and the 
properties of the characteristic directions discussed 
previously are of course retained. Furthermore the 
formulation presented in equation (32) is uniformly 
valid for subsonic, transonic, and supersonic flow 
conditions. 



This is a general surrogate second-order system for 
equation (29). Its properties depend of course on the 
specific choices made for the 4x4 matrices Af and N. 

In particular, if we choose M—NAB~^ and A as 
some nonsingular constant matrix, such as N=I, the 
identity matrix, then equation (30) reduces to the 
particularly simple form 



or more precisely 


Summation 

We have demonstrated the capability of the 
surrogate-equation technique to generate second- 
order partial differential systems suitable for use in 
transonic flow computations. The fully conservative 
algorithm developed here to solve the surrogate 
second-order system for the transonic small- 
disturbance equations appears to be superior, when 
measured in terms of computational efficiency, to 
the standard Murman and Cole approach to solving ^ 
the perturbation-potential equation. 

The surrogate-equation technique shows promise 
toward providing a method for iteratively solving the 6. 

full steady Euler equations in the transonic flow 
regime. Indeed, the results obtained thus far indicate 
that further research and algorithm development 
aimed at producing such a method are merited, and 
such work is presently in progress. 

Finally, although the present discussion has been 
limited to two dimensions, the eventual extension of 
the surrogate-equation technique to three dimensions 
appears to be entirely feasible. 
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Summary of Results 

We have been concerned during the course of this 
research with the development of a method to obtain 
approximate numerical solutions to fairly general 
classes of inviscid flow problems. We have chosen, 
on the one hand, to avoid the usual potential and 
stream function approaches with their inherent 
limitations and, on the other hand, to try to devise a 
method capable of greater computational efficiency 
than the time-asymptotic solution of the unsteady 
Euler equations. 

The search for such a method has led to what we 
refer to as the surrogate-equation technique. This 
technique provides a means for obtaining, from a 
given first-order partial differential system, a 
replacement second-order system whose solution set 
contains the sought-after solution to the first-order 
system. The possibility of other, undesirable, 
solutions being admitted is eliminated by using the 
original first-order system to supply the additional 
boundary conditions required by the surrogate 
second-order system. Given the surrogate second- 
order system, it is a relatively straightforward matter 
to construct for it stable, fully conservative, type- 
dependent finite difference equations that can then 
be solved by using, for example, a successive-line- 
relaxation iterative procedure. 

The essential aspects of the surrogate-equation 
technique were illustrated earlier, where it was 
applied to the Cauchy-Riemann equations. There 
numerical solutions of the surrogate second-order 
system were compared with those of the standard 
velocity potential formulation. 

Problems of greater complexity, Euler problems, 
were chosen so as to more thoroughly display the 
utility and generality of the surrogate-equation 
technique. The full, steady Euler equations were used 
to model the rotational, inviscid subcritical flow 
through a two-dimensional bend. The relaxation 
solutions of the surrogate .second-order system were 
in this case contrasted with analytic, incompressible 
flow solutions obtained by Shercliff. 

The question of the applicability of the surrogate- 
equation technique to transonic flow computations 
was addressed next. A surrogate second-order system 
was constructed for the transonic small-disturbance 
equations. A solution algorithm based on the 
surrogate-equation formulation was then compared 


with the standard Murman and Cole type of 
algorithm for the perturbation-potential equation. 
The surrogate-equation algorithm produced good 
results and surpassed the computational efficiency of 
the perturbation-potential algorithm by several 
multiples. An indication was also given as to how the 
transonic results might be extended to the full, steady 
Euler equations. 

The main conclusions are as follows: 

1. It is possible to obtain an approximate 
numerical solution to a system of first-order partial 
differential equations by solving a problem 
consisting of a surrogate second-order system 
together with the original boundary conditions and 
supplementary relations obtained from the first- 
order system. 

2. The surrogate-equation technique can be 
applied successfully to the solution of inviscid flow 
problems across the entire spectrum of subsonic, 
transonic, and supersonic conditions. 

3. The surrogate-equation technique provides a 
means for formulating problems involving the first- 
order equations describing inviscid flow in such a 
way as to allow the use of fully conservative, type- 
dependent differencing and iterative solution 
procedures. 

4. The surrogate-equation technique provides a 
means for solving inviscid flow problems without 
resort to assuming the existence of either a velocity 
potential function or a stream function. Hence the 
surrogate-equation approach is more general than 
either of these other approaches and does not share 
their limitations. In particular, the surrogate- 
equation technique can be employed to obtain 
solutions to flow problems where any combination of 
rotationality, transonic conditions, or three 
dimensionality is present. 

5. An application of the surrogate-equation 
technique to the transonic small-disturbance 
equations results in an algorithm that, on the basis of 
the computational experimentation reported herein, 
appears to have a computational efficiency several 
times greater than that of the standard perturbation- 
potential algorithm. 


Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, November 19, 1980 



Appendix A 

Jacobian Matrices and Some of Their Properties 

Recall that if the vector of unknowns w is defined such that w = \p,pu,pv, E] the Jacobian matrices A and 
B are defined as >1 = dfidw, where f and ^ are as given in equation (10) of the main text. The matrices A and B 
can be written explicitly as 
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uv 

— V 

— u 

0 

3-7 , 1 -7 2 

( 7 - 1 )m 

(7-3)t; 

1 -7 

+(1 —y)v(u^ + v^) 
P 

( 7 — l)uu 

P 2 

— yv 


To examine some of the properties of A and B, we first apply the similarity transformations S and Sj , as 
given by Turkel (ref. 50; see also refs. 51 through 53) to create the matrices A j and Bi , which are similar to A 
and B, respectively. 
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Ai =SiSAS-^Si~' 
Bj =SiSBS-^S,-' 


These transformations yield 
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Hence it is clear that A is singular for «= ±c and 
u = 0and£is singular for v- ±cand v = 0. Also the 
eigenvalues of BA are clearly 
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Appendix B 

Simplified Expression of Euler Equations 


Given the Euler equations in the form 
fx + {TF)y= 0 , where T=BA~^, we can carry out 
the indicated differentiation to arrive at 

fx + Tfy + Tyf= 0 

Since g=Tf, we can write (g — Tf)y = 0, which can 
then be reexpressed as 


Sy-Tfy-Tyf =0 


But it follows immediately that gy — Tfy =0 and 
hence that Tyf= 0 . Thus the Euler equations can be 
written in simplified form as fx-^Tfy=0, which 
appears as equation (11) in the main text. 

Similarly Txf=0 and consequently 

Txfy = Tyfx 

This latter observation allows reexpression of the 
surrogate second-order system for the Euler 
equations 


dx^ dy^ 


(TTy-Tx)^ 


^=0 


as 


a2 ;)2 a 

-^-T^-^-iTTy + TyT) — 

■ dx^ dy^ . 


^=0 


This equation can then be rewritten in the form 


-dx^ 



which appears as equation (14) in the main text. 
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